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Abstract. We demonstrate that the invaded cluster algorithm, recently introduced 
by Machta et al, is a fast and reliable tool for determining the critical temperature and 
the magnetic critical exponent of periodic and aperiodic ferromagnetic Ising models 
in two dimensions. The algorithm is shown to reproduce the known values of the 
C^ , critical temperature on various periodic and quasiperiodic graphs with an accuracy 

S ' of more than three significant digits, but only modest computational effort. On two 

' ■ quasiperiodic graphs which were not investigated in this respect before, the twelvefold 

(-H , symmetric square-triangle tiling and the tenfold symmetric Tiibingen triangle tiling, 

Q ' we determine the critical temperature. Furthermore, a generalization of the algorithm 

, ^, . to non-identical coupling strengths is presented and applied to a class of Ising models 

on the Labyrinth tiling. For generic cases in which the heuristic Harris-Luck criterion 
Cn ■ predicts deviations from the Onsager universality class, we find a magnetic critical 

i*Zi , exponent different from the Onsager value. But also notable exceptions to the criterion 

fT^ ' are found which consist not only of the exactly solvable cases, in agreement with a 

recent exact result, but also of the self-dual ones and maybe more. 
C3 ■ 
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'^ ! 1. Introduction 

o 

O ' Explicit exact solutions for the thermodynamic behaviour of (ferromagnetic) Ising 

models are available, or expected to exist, only in one dimension and for a few models in 
^ ■ two dimensions. For more complex cases, and in higher dimensions, one has to rely on 

j^ ■ approximate techniques or simulations. As one is usually interested in the expectation 

value of certain observables such as the magnetization and the susceptibility at a given 
temperature, Monte Carlo algorithms sampling the canonical ensemble are often the 
numerical methods of choice. 

The limiting factor in investigating the critical behaviour of spin systems with 
algorithms of this kind is the critical slowing at the phase transition point. The 
time needed to sample statistically independent system configurations diverges with the 
system size according to a power law. The algorithms being most efficient in this respect 
are the so-called cluster algorithms, the first of which was introduced by Swendsen and 
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Wang in 1987 IQ. A single-cluster version, which shows less critical slowing in three 



and more dimensions, was proposed by Wolff two years later p6[. Also series expansions 



are suitable for studying the critical behaviour of periodic models, but they do not lead 



to satisfactory results for aperiodic ones pO 



Estimating the critical temperature of the infinite system with Monte Carlo 
algorithms is usually done using the Binder cumulant method [0. (E.g. this has been 
used in the first paper on the universality class for quasiperiodic graphs f^ which will 
play an important role also in this paper.) For a number of finite system sizes and several 
temperature values taken from some interval containing the critical temperature, the 
fourth cumulant of the magnetization is measured. This cumulant is zero for infinite 
temperature, | for zero temperature and depends on the system size for all other 
temperature values; only at the critical temperature the curves for different system sizes 
intersect. This admits rather reliable determination of the location of the critical point. 
However, because the system has to be simulated at several temperature values, this 
procedure is very time-consuming, and a rough knowledge of the location of the critical 
temperature is required beforehand. So, independent methods are certainly useful. 

In 1995 Machta et al p4| developed a self-organized version of the Swendsen-Wang 



algorithm which they dubbed invaded cluster (IC) algorithm. It is not only able to locate 
the critical temperature without prior knowledge and without any temperature sweeps 
being necessary, but also seems to show yet less critical slowing than the Wolff algorithm 
p6| . The algorithm was shown to reproduce three significant digits of the known values 
of the critical temperature on the square and simple cubic lattice Ising models with 
modest computational effort [Bl, ES . The ensemble that is sampled, however, is not the 



canonical one, and not much is known about it rigorously. In the thermodynamic limit, 
it is expected to be equivalent to the latter, but about the finite-size scaling ^ nothing 
helpful is known. Also, the thermal exponent and equivalent quantities, such as z/ and 
a, are not easily measured with the IC algorithm ||2^, [l^ . 

This lack of theoretical knowledge is responsible for the fact that so far the best 
accuracy in measuring the critical temperature is still an order of magnitude lower than 
with other cluster algorithms [0. Nevertheless, at this accuracy, the IC algorithm is 
substantially faster. This is mainly due to the fact that the temperature need not be 
varied. In addition, the autocorrelation times are by a factor of 4 or more smaller than 
those of the Wolff algorithm p5| . (If one is also interested in observables on sub-systems 
this is not true, but the autocorrelation times are still comparable to the ones of the 



Wolff algorithm, see |26].) As a rule of thumb, one can expect the IC algorithm to be at 
least a factor of 5 faster. Therefore it is a good choice for obtaining a good estimate of 
the critical temperature with modest computational effort. In addition, one of the two 
independent critical exponents, the magnetic exponent yh, can be measured easily. 

Those quantities are of particular interest for Ising models on quasiperiodic tilings, 
a class of graphs being used as models for the structure of quasicrystals, compare 0, ^ 
and references therein. As these graphs are very homogeneous, the influence of the 
quasiperiodicity on the critical temperature is expected to be small and its rough location 
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should be governed by the first mean coordination number, see 0] and references therein. 
The critical exponents are expected to be the same as their periodic counterparts 
unless the degree of disorder is bigger than a 'critical' value. This is the result of 
the heuristic criterion by Luck |^ who generalized the well-known Harris criterion for 
random disorder |T6[ to aperiodic structures. This criterion, now called the Harris-Luck 
criterion, was recently proved for certain kinds of one-dimensional disorder ||2y, T7| and 
corroborated by approximate methods for two-dimensional substitution systems [|l^, [T^ . 
But in two or more dimensions, the validity has not yet been systematically confirmed 
by simulations. 

The IC algorithm has so far only been tested on the square and the simple cubic 
lattice - both with success 123, |25|, O, |i 



||. Before its results for more complex situations 
can be trusted, however, it should be tested on other graphs for which the critical 
temperature is known. This will be done in section ^ for a number of periodic and 
quasiperiodic graphs. Then, as a first application, the critical temperature of the Ising 
model on two quasiperiodic tilings, the twelvefold symmetric square-triangle tiling [^ 
and the tenfold symmetric Tiibingen triangle tiling , will be determined. In section Q, 
we will present a generalization of the IC algorithm to models with arbitrary coupling 
strengths. Results of simulations on different realizations of the Labyrinth tiling [^ 



with this generalized algorithm will be described in section ^. Section ^ will give a 
summary and discussion. 

2. The critical temperature of models with identical couplings 



We consider the field-free Ising model defined by a spin ai = ±1 on each vertex i of some 
graph and a ferromagnetic bond {i,j) of coupling strength Jij > between each pair 
of neighbouring spins. The internal energy is ?^(cr) = —J2{i j) JijO'iO'j where the sum is 
over all bonds. In this section, we restrict ourselves to identical couplings Jij = 1. 

For the Ising model with this restriction, the critical temperature is known exactly 



on a number of periodic graphs, see |^ for a survey. In the case of two quasiperiodic 
tilings, the Penrose tiling p8[ and the octagonal Ammann-Beenker tiling [|l], ^, there 
exist Monte Carlo estimates ^1, 33, 21] and high-precision numerical values from a 
recent analysis of large periodic approximants using Kac-Ward determinants pTl. With 



those values, also the critical temperature on the corresponding dual tilings is known 



through the exact relationship |35 



sinh(2/3e) sinh(2/?;) 



(1) 



where /3c and (3* are the inverse critical temperatures on a graph and its dual. All this 
offers a number of graphs on which the IC algorithm can be tested thoroughly. 

With the IC algorithm, a given configuration of a (finite-sized) system is updated in 
two steps. First, bonds between aligned spins {satisfied bonds) are selected in random 
order until one connected component (cluster) wraps around the system in at least one 
dimension (periodic boundary conditions are assumed) or all satisfied bonds have been 
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Figure 1. Parts of the periodic graplis considered here: (a) square lattice, (b) 
triangular lattice, (c) hexagon packing, (d) Kagonic graph, and (e) diced graph. 





Figure 2. Parts of the quasiperiodic tilings considered here: (a) Penrose, (b) 
Animann-Beenker, (c) Tiibingen triangle, and (d) twelvefold square-triangle tiling. 
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selected. Then, each cluster (including each isolated spin) is flipped with probability 
i. The fraction of selected to satisfled bonds, /, gives an estimate /3est of the inverse 



critical temperature via pi 



/ = l-e-2/3-'. (2) 

This is reasonable as the expectation value of / in the canonical ensemble is just the 
probability p = 1 — e'"^^ to select a bond in the Swendsen-Wang algorithm. 

For each periodic graph depicted in flgure |I|, we used the IC algorithm to simulate 
the Ising model at 13 different linear system sizes L between 16 and 288. For the 
quasiperiodic tilings, shown in flgure |^, we had to restrict ourselves to the periodic 
approximants available. With each size, 10 bunches of 5,000 update steps each were 
taken. This took approximately 2.3 ■ 10~^ seconds per bond on a 266 MHz Pentium 
II processor which amounts to roughly one day for all sizes of one graph. The flrst 
bunch was used for equilibration and then discarded. For each one of the remaining, 
the averages of the quantities of interest were computed and those numbers treated as 
an independent sample. 

In contrast to the previous IC studies fl^ , ^, ^ |26[ , in which the average of / was 
determined and /3est inferred in the end using (H), we computed /3est in each step. This 
leads to the same results after flnite-size scaling, as shown in [^. Figure |^ shows the 
results for the estimates of the inverse critical temperature on the hexagon packing as a 
typical example. The picture is qualitatively the same for all graphs considered, details 
can be found in |^. Due to certain bottlenecks occuring in the cluster growth process 
p5| , the mean value of /3est is always greater than the median. But both are expected 
to approach the same value /3c in the thermodynamic limit L —>■ oo. From the data 
it is obvious that both curves are inclined towards each other. This gives rise to the 
following recipe for estimating (3^ and its error. We compute linear flts to both the mean 
and median values independently using the data points corresponding to some of the 
largest system sizes. Then we take the mean of both /5c-axis intercepts as an estimate 
for /?c and half their separation as its error. 

There is no a priori justiflcation for this recipe besides the numerical evidence that 
/3c lies in between both points. The results shown in table |l| agree with the known values 
to within more than three signiflcant digits. This indicates that the way /3c is estimated 
is reasonable and that we considerably overestimate the error. For the Penrose and the 
Ammann-Beenker tiling, the agreement of our data with the high-precision values [^] is 



even better than the one of the other Monte Carlo values [^, ^, ^ . This is especially 
remarkable if one takes into account the small computational effort. 

One can conclude from this that the values of the critical temperature for the 
twelvefold symmetric square-triangle tiling [^ and the tenfold symmetric Tiibingen 
triangle tiling ^ given in table [^ should have a similar accuracy. Also, they are 
in the correct rough range that one expects from the comparison of their flrst mean 
coordination number with the one of the square and triangular lattice, see table |l[ 

The autocorrelation times found were all compatible with the ones given in table II 
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Figure 3. Finite-size scaling of the IC estimates of the inverse critical temperature 
on the hexagon packing. 



Table 1. IC estimates of the inverse critical temperature for various graphs. Note 
that the errors given are very conservative upper bounds. Exact and other numerical 
values are given for comparison where available. The exact values are taken from |35|] , 
the best numerical ones from |3^. The values for the dual tilings marked with * were 
inferred using the exact relation (yj). 



graph 




/3c 






(mean) 

coordination 

number 






exact/best 


IC estimate 


other values 






square lattice 


0.4406868 . . . 


0.4405(4) 






4 




triangular lattice 


0.2746531 . . . 


0.2747(5) 






6 




hexagon packing 


0.6584789 . . . 


0.6584(11) 






3 




Kagome graph 


0.4665660 . . . 


0.4665(8) 






4 




diced graph 


0.4157215... 


0.4157(7) 






4 




Penrose tiling 


0.417046(1) 


0.4170(8) 


0.4181(7) 1 
0.4165(9) 


27 
33 


1 4 
1 




dual Penrose tiling 


0.465145(1)* 


0.4652(8) 






4 




Ammann-Beenker tiling 


0.41887800(1) 


0.4191(7) 


0.4186(7) 1 


21 


1 4 




dual Ammann-Beenker 


0.46318974(1)* 


0.4634(10) 






4 




Tiibingen triangle tiling 




0.2565(6) 






6 




square-triangle tiling 




0.3364(5) 






5.072...- 12- 
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of [Q, especially they were always less than unity. Therefore it seems safe to conclude 
that the IC algorithm is not significantly slower on quasiperiodic graphs than on periodic 
ones. In particular, it should still be much faster in finding the critical temperature (to 
the above mentioned accuracy) than other Monte Carlo algorithms. 

3. Generalization of the IC algorithm to non-identical coupling strengths 

The IC algorithm was so far only published for identical coupling strengths Jij = J, 
but a generalization to arbitrary ones is rather straightforward |^. To this end, let 



us review how the Swendsen-Wang (SW) and IC algorithms are connected. In the SW 
algorithm, for each satisfied bond {i,j) a random number < tij < 1 is drawn from a 
uniform distribution and the bond is selected if the condition 

tij < Pi, := 1 - e-^^^^^ (3) 

is fulfilled. In the IC algorithm, on the contrary, bonds are selected in random order 



until one cluster wraps around the system (or a different condition is met, see ||2^, |T2[)- 
After solving (0) for /?, 

/3<--^log(l-t,,)=:A„ (4) 

we see how the algorithm has to be generalized. We again draw a random number tij 
for each bond and calculate the corresponding f3ij from (^). Then we sort the bonds 
ascendingly with regard to Pij and select the satisfied ones in this order until one cluster 
wraps around the system (or all satisfied bonds have been selected). If bond {k,l) was 
selected last, Pki gives an estimate /?est for the inverse critical temperature. 

The random numbers ty play a twofold role here. They determine the order in 
which the bonds are selected, and the one corresponding to the last bond selected gives 
the numerical value of /3est- So it is not necessary to use independent distributions for 
the tij as long as the marginal distributions are the same for all of them. Instead, one 
can use a random permutation of equidistant points 'approximating' the interval [0, 1[ 
for the tij. With such a choice, the original algorithm is recovered for identical couplings 
(in which case the sorting is trivial). 

Now, there is a crucial difference between the original and the generalized version 
with regard to the computational complexity. In the case of identical couplings, only a 
permutation of the A^ bonds has to be created. This can be done with a computational 
effort of order O(A^). The sorting of the arbitrary numbers (3ij, in contrast, requires an 
effort of O(A^logA^) which makes the algorithm considerably slower. If, however, only 
a small (constant) number k of different coupling strengths is present, the sorting can 
be done with effort 0{kN) = 0{N) by the following steps. (Let us assume all bonds to 
be numbered sequentially 1, . . . , A^ here.) 

1. Create an auxiliary permutation n of the numbers {1, . . . , A^}. 
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2. For each bond, define a value of tlie inverse temperature, 

/?.(.) :=-^ log (l-^), ^ = l,...,N. 

3. Create tlie final permutation 11 of {1, ... , N} by repeating the following steps for 
i = 1, . . . ,N after initializing a list {kt)i<t<k with the indices of the first bonds of 
type t appearing in vr. 

(a) Determine the value of t for which Pnikt) becomes minimal. 

(b) Set n(z) := 7r{kt). 

(c) Replace kf by the next bond of type t appearing in vr. 

Steps 1 and 2 are possible with effort O(A^). In step 3, each of the k counters kt takes 
on all the values {1, . . . , A^} in the worst case. This step can therefore also be performed 
with effort 0{kN) = 0{N). Overall, the sorting procedure takes linear computational 
effort, only with a higher constant than in the original algorithm. 

The generalized IC algorithm will be applied to a class of Ising models with non- 
identical coupling strengths in the following section. 

4. The Labyrinth tiling as an example to study the relevance of disorder 

4.I. The Ising model on the Labyrinth tiling 



The Labyrinth tiling |3^ is constructed from a substitution p on the two-letter alphabet 



A = {a, b}. Here we will consider substitutions of the type 



a ^ b 

b -^ baH ' 



P • L 7,^fc 



with k>l. The substitution is iterated infinitely, starting from the word b (or any other 
finite word), i.e. in each step it is applied to each letter of the word. This yields a (semi-) 
infinite limiting word w. For A; = 1, we have the silver mean chain treated in [Q. Each 
letter is assigned an interval of a specific length. In this representation, the Cartesian 
product of w with itself is taken giving one quadrant of a (distorted) square lattice. Of 
this, one takes every other vertex, starting from the origin, and connects each one with 
its nearest neighbours over the diagonals of the underlying cells. As there are four types 
of those cells according to the letters x,y E A on the vertical and horizontal copy of w, 
and each bond can be either 'raising' or 'lowering', we can distinguish 8 different types 
of bonds in the corresponding Ising model. We will denote the corresponding coupling 
strengths by J^y for the raising and Jxy for the lowering bonds. Part of the Labyrinth 
tiling corresponding to fc = 1 is shown in figure ^. 

On a 4-dimensional submanifold of the 8-dimensional coupling space, the Ising 
model is exactly solvable and shows critical behaviour of Onsager type regardless of the 
underlying word w [If. To see what the Harris-Luck criterion [|l^, |2^ predicts, we have 



to look at the fiuctuation exponent a;, which describes how the deviations of the mean 
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Figure 4. Part of the Labyrinth thing corresponding to the substitution with k = \. 



coupling strength in a finite patch, compared to the infinite volume mean, scale with the 



size of the patch. If uj is greater than the 'critical' value uOc = 1 — ^|— j; ||2^, one expects 
the disorder to be relevant, i.e. critical exponents deviating from the Onsager values. 
Here, v is the correlation length exponent of the periodic system (i/ = 1 in our case) 
and the disorder affects all d^ = 2 dimensions, thus 0Jc = \. For oj < Uc, the disorder is 
expected to be irrelevant. In the case of marginal disorder, u = Uc, the criterion does 
not make any specific predictions. 

For substitution systems, the fiuctuation exponent u: can be extracted from the 
bond substitution, which is in our case induced by the letter substitution p. It is given by 
uj = log |A2|/logAi, where Ai > IA2I are the two largest eigenvalues of the corresponding 
substitution matrix whose entries count how many bonds of one type (determined by 
the row of the matrix) are in the substitute for another type (the column). This 
matrix turns out to be described by the tensor product of the substitution matrix 



'Ok' 

,12, 



with a{Mp 



for p with itself, its spectrum being a{Mp ® Mp) = {A/x : A,/i G a{Mp)} 
= {1 ± vT+fc}. Thus, we find the fiuctuation exponent 

UJ = , (5) 

2iog(i + ym) ^^ 

Substitutions with A; < 3 correspond to irrelevant disorder. A; = 3 is the marginal case, 

and for A; > 3 the disorder should be relevant. These predictions shall be tested shortly. 

As an alternative model, let us also look at a random case. Instead of using a 

substitution for creating the word w, one can also use a random word of which each letter 
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Table 2. Choices for the couphng strengths used in the simulations for the Labyrinth 
tiUng. The ones for the raising and lowering bonds are denoted by Jxy and J^y, 
respectively, where x and y are the corresponding letters in the horizontal and vertical 
copy of the word w. 



choice 


'^aa 


Jab 


Jba 


Jbb 


"^aa 


Jab 


Jba 


Jbb 


#1 


1.61268 


1 


0.63693035 


0.35675743 


0.55876724 


1 


1.4707784 


2.1084272 


#2 


1.3403778 


1 


0.77308147 


0.55051591 


0.71949498 


1 


1.2644651 


1.6289055 


#3 


0.4 


0.8 


1 


2.8 


1.9809928 


1.2285752 


1 


0.19281531 


#4 


2.5 


1 


0.5 


0.75 


0.25159591 


1 


1.7343053 


1.2964035 


#5 


1.2 


1 


0.7 


2 


1.4 


0.9 


1.1 


1.8 


#6 


2.5 


1 


0.4 


0.8 


2.8 


0.9 


1.3 


0.7 


#7 


0.4 


0.8 


1 


2.8 


1.3 


2.3 


0.9 


2.4 


#8 


2 


1.2 


1.2 


0.5 


2 


1.2 


1.2 


0.5 



is independently chosen to be a or 6 with probabihty pa and (1 — Pa), respectively. Due 
to the way the Labyrinth is constructed, the resulting bond distribution has correlations 
within each row and column. Accordingly, this kind of disorder has a higher fluctuation 



exponent than uncorrelated disorder, u 
respect to the Harris-Luck criterion. 



compared to uj 



|, and is relevant with 



4-2. Results 

First, we successfully tested the implementation of the generalized algorithm on graphs 
with identical couplings and a few of the exactly solvable cases of the Labyrinth tiling. 
Then, we applied the IC algorithm to the Ising model on the Labyrinth 
corresponding to the substitutions k = 2, 3, 4 and on the random version with pa = 0.4, 
using the choices for the coupling strengths given in table 0. The flrst two choices 
(#1, #2) are on the exactly solvable submanifold, i.e. we chose three of the coupling 
strengths {Jaa, Jab, Jba) arbitrarily and the inverse critical temperature to be the one of 
the square lattice Ising model, /3c = arsinh(l)/2 ~ 0.44069, and determined the other 
flve couplings by numerically solving the flve equations H 



and 



sinh(2/3c Jxy) sinh.{2(3c Jxy) = 1, x,y e A, 



Sinh(2/5c Jaa) sinh(2/?cJbb) V COsh(/?c[- Jaa + Jba + Jbb + Jab]) 



(6) 



Sinh(2/5c Jafe) Sinh(2/?ejfea) / COsh{/3,[-Jaa + Jba + Jbb + Jab]) 
cosh (/3c [Jag - Jba + Jbb + Jab]) COsh{Pc[Jaa + Jba - Jbb + Jgb]) 
COsh(/3c[Jaa " Jba + Jbb + Jab]) COsh(/3c[Jaa + Jba " Jbb + Jab]) 
cosh (/3c [Jaa + Jba + Jbb " Jab]) -. 



(7) 



cosh (/3c [Jaa + Jba + Jbb " Jab]) 

The flrst four (^ are just the duahty conditions on the lowering and raising bonds of 
the same type. The flfth one (0) involves all the remaining four coupling strengths. 
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0.345... 
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1.876 
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1.874 


1.875 


1.873 
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0.590... 


1.872 


1.873 


1.877 


1.875 


random 


3 

4 


1.874 


1.878 


1.877 


1.862 
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Table 3. IC estimates of the magnetic critical exponent yh on the Labyrinth tiling 
for three different substitutions and the random case, each one with the eight choices 
of the couphng strengths given in table 0. The value marked with * is significantly 
different from the Onsager value y^ = 1.875. Due to fluctuations as visible in part (b) 
of figure o, the total error for the values marked with * or f is presumably quite large 
(« 0.03). In the random case, the flucuations for choices #5 to #8 were too large to 
even confirm a power-law behaviour. The fiuctuation exponent uj for the substitutions 
was computed using (^. 

choice of couplings 

k u; #1 #2 #3 #4 #5 #6 #7 #8 

1.875 1.871 1.859 1.872 
1.871 1.847t 1.732t 1.880 
1.866 1.839t 1.453* 1.898t 
fluctuations too large 



The next two choices (#3, #4) are self-dual but not exactly solvable, i.e. they fulfill the 
duality conditions (|^) with respect to our chosen temperature, but the fifth condition 
(0) is significantly violated. The last four choices (#5 to #8) violate all five conditions. 
Choice 7^8 is isotropic, i.e. bonds are only distinguished with regard to their length. 

We estimated the magnetic critical exponent yh = 2 — [3 /v as described in [Q 
by plotting the logarithm of the expectation value of the mass M of the cluster that 
wrapped around the system against the logarithm of the linear system size L and then 
taking a linear fit. As the expectation values of M and the magnetization of the system 
are the same, the usual scaling relation for the magnetization becomes M oc L^^ . Two 
examples are shown in figure ^ the results are summarized in table ^. For the exactly 
solvable choices (#1, #2), we always found the exponent to be of Onsager type and the 
critical temperature to agree with the exact solution, i.e. with the choice made above. 
For the self-dual choices (#3, #4), again, the Onsager exponent was found in all cases. 
The critical temperature, however, was not the chosen one for even values of k (2, 4), 
but agreed with it for odd ones (1, 3, 5) and probably also in the random version. (For 
the latter, the fiuctuations in /3est were too large to state this unambiguously.) This is 
not yet understood. 

For the arbitrary choices (#5 to #8), the critical exponent was of the Onsager 
class for k = 2. For A; > 3, we still found power-law behaviour for all choices, but a 
significantly different exponent for choice #7 with k = 4. For the values in table ^ 
marked with f, no unique conclusion could be drawn, although the linear fits indicate 
deviations. With the clear deviation for choice #7, k = 4, and the predicition of the 
Harris-Luck criterion in mind, also choice #7 with k = 3 might be interpreted as a 
deviation from the Onsager value. Choice #5 and choice #8 with k = 2 are definitely 
compatible with the Onsager exponent, but the deviations might just be too small to 
be detected. In the random case, the fiuctuations in M were too large to even confirm 
a power-law behaviour. This is probably due to the random fiuctuations of the bond 
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Figure 5. Determination of tlie magnetic exponent yh on tlie Labyrintli tiling. Tlie 
mass M of the cluster wrapping around the system is plotted versus the linear system 
size L. (a) is a typical example of Onsager behaviour (coupling choice ^4 with 
k — 4) and (b) is typical of non-Onsager behaviour (choice ^7 with fc = 4). The 
solid lines have the Onsager slope yh — 1.875, the dashed line is a linear fit yielding 
yfi ~ 1.453. The fluctuations in (b) are due to variations in the bond frequencies in 
the approximants used. 
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frequencies in the approximants used. For choices 7^1 to #4, however, the picture was 
quahtatively the same as in part (a) of figure |^. 

5. Summary and discussion 

We tested the IC algorithm for the Ising model on various periodic and aperiodic planar 
graphs. A procedure was described that makes use of special properties of the algorithm 
in two dimensions and allows to determine the critical temperature with an accuracy of 
more than three significant digits. The computational effort for this is small compared 
to other Monte Carlo methods. Then we estimated the critical temperature on the 
twelvefold symmetric square-triangle tiling and the tenfold symmetric Tiibingen triangle 
tiling, two cases for which no values were known before. The values found are located 
in the rough range that one expects from the first mean coordination number of the 
tilings. 

In the second part of the paper, we presented a generalized version of the algorithm 
applicable to models with non-identical coupling strengths. We applied it to the Ising 
model on the Labyrinth tiling for three different underlying substitutions, corresponding 
to irrelevant, marginal, and relevant disorder according to the Harris-Luck criterion, and 
a random case with relevant disorder. Each one was simulated for a few typical choices 
of the coupling strengths on the exactly solvable submanifold, on the self-dual but not 
exactly solvable submanifold, and away from both. The magnetic critical exponent 
Hh = 2 — (3/v was determined in all cases. The values found were compatible with the 
Onsager value for all self-dual subcases including the exactly solvable ones. For the 
latter this was known to be true exactly |^. But it seems to extend to all self-dual 
cases even when the Harris-Luck criterion predicts deviations. For some (but not all) 
other choices of the coupling strengths, when the disorder was marginal or relevant 
according to the criterion, power-law behaviour with different exponents was observed 
in the substitution systems. The criterion can therefore be expected to be generically 
correct, but it does not exclude exceptions on lower- dimensional coupling manifolds. 
The sharp contrast seen in the random version between self-dual and not self-dual cases 
with respect to the absence, respectively presence, of large fluctuations of the mass of 
the cluster that wraps around the system might indicate that exceptions to the criterion 
are restricted to the self-dual submanifold. But further investigations are necessary to 
confirm this observation. 

The IC algorithm has proved a good tool for efficiently obtaining quite accurate 
estimates of the critical temperature for periodic and aperiodic planar graphs. Although 
only one of the two independent critical exponents, the magnetic one, can be measured, 
this gives at least the possibility of detecting deviations from the Onsager universality 
class. Fortunately, it is a more sensitive quantity for this purpose than the specific heat 
exponent a, which is expected to remain zero when the disorder is increased [O. 
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